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Abstract 

In quantum theory, the "spinless Salpeter equation," the relativistic generalization of the 
nonrelativistic Schrodinger equation, is used to describe both bound states of scalar particles 
and the spin-averaged spectra of bound states of fermions. A numerical procedure solves the 
spinless Salpeter equation by approximating this eigenvalue equation by a matrix eigenvalue 
problem with explicitly known matrices. 



PACS numbers: 03.65. Ge, 03.65.Pm, ll.10.St 



* E-mail: wolfgang.lucha@oeaw.ac.at 
f E-mail: franz.schoeberl@univie.ac.at 



1 



1 Introduction 

The appropriate framework for the description of bound states within a relativistic quantum 
field theory is, beyond doubt, the famous Bethe-Salpeter equation [jij. Its actual application, 
however, faces several well-known problems of both conceptual as well as practical nature. 
For this reason, one usually considers some (three-dimensional) reduction of this formalism. 
In particular, assuming a static interaction between the bound-state constituents, neglecting 
all spin degrees of freedom, and considering exclusively positive-energy solutions, one ends up 
with what is called the "spinless Salpeter equation." For some brief accounts of the reduction 
of the Bethe-Salpeter equation to the spinless Salpeter equation, see, for instance, Refs. |2|, |3|. 
This eigenvalue equation generalizes the Schrodinger equation of the nonrelativistic quantum 
theory towards relativistic kinematics. It thus accomplishes the description of bound states of 
spinless constituents (scalar bosons) as well as of the spin-averaged spectra of bound states of 
fermions, like, e. g., in elementary particle physics the description of hadrons as bound states 
of quarks Q within the (intuitive) framework of potential models for the strong interactions. 

Here we present a very efficient method, based on Ref. ||, for solving the spinless Salpeter 
equation in configuration space, where the interaction potentials are usually formulated. The 
efficiency of this method stems from the fact that it does not require numerical integrations; 
rather, the aim of this method is to employ analytical results whereever these are available. 

The outline of this paper is as follows. In Sect. |2|, the configuration-space representation of 
the spinless Salpeter equation is briefly recalled. The resulting integro-differential equation is 
then approximated, in Sect. ^, by a (finite-dimensional) matrix eigenvalue problem. The only 
nontrivial ingredient of this matrix equation, viz., the piece emerging from the kinetic energy, 
is calculated in Sect. ||. This procedure for the determination of both energy eigenvalues and 
corresponding wave functions of bound states described by a spinless Salpeter equation with 
an, in principle, arbitrary interaction potential is eventually condensed, in Sect. g[ to a simple 
Mathematica routine. (This routine may be obtained, of course, by contacting the authors.) 

2 The Spinless Salpeter Equation in Configuration Space 

For notational and conceptual simplicity, we consider only the case of equal masses m of the 
bound-state constituents; the generalization to the case of unequal masses is straightforward. 

The semirelativistic Hamiltonian H governing the dynamics of two relativistically moving 
spinless particles of equal masses m, which interact via some arbitrary coordinate-dependent 
static potential V = V(x) reads, in the center-of-momentum frame of these two particles, 

H = 2T + V, (1) 

where T denotes the "square-root" operator of the relativistic expression for the free (kinetic) 
energy of some particle of mass m and momentum p, 

T = T(p) = ^/p 2 + m 2 . (2) 

The spinless Salpeter equation is nothing else but the eigenvalue equation for the operator H, 

H\ Xk )=E k \xk) , k = 0,1,2,... , (3) 

for the complete set of Hilbert-space eigenvectors \xk) with corresponding energy eigenvalues 

T-i _ (Xk\H\xk) , n 1 o 

E k = —, — i — r~ > fc = 0,1,2,.... 

\Xk\Xk) 

However, in contrast to, e. g., the (nonrelativistic) Schrodinger equation, the semirelativistic 
Hamiltonian H is a nonlocal operator, i. e., either the relativistic kinetic-energy operator T in 
configuration space or, in general, the interaction-potential operator V in momentum space 
is nonlocal. 
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In configuration space, the action of the kinetic-energy operator (J2j) on some element if) of 
L2(i? 3 ), the Hilbert space of square- integrable functions on the three-dimensional Euclidian 
space R 3 , is defined by 



(TVO(x) 



(27T) 



— rg y d 3 p y d 3 y \jv 2 + m 2 exp [i p • (x 



y)]V>(y) 



With this definition, an integral representation of the spinless Salpeter equation can be found 
H |[|. We focus our interest, for the moment, to spherically symmetric interaction potentials 
V(x) = V(r), depending only on the radial coordinate r = |x|. For a generic eigenfunction x 
of H describing a state of definite orbital angular momentum £, we introduce a reduced radial 
wave function w(r) by factorizing x( x ) according to 

X (x)=r i w(r)y em (n) , 

where 3^ m (f2) denote the (orthonormalized) spherical harmonics of angular momentum t and 
its projection m, depending on the solid angle f2. For a nonvanishing mass of the bound-state 
constituents, that is, m > 0, upon rescaling the radial variable r according to x := mr and 
introducing the scaled counterpart w{x) of the reduced radial wave function w(r) by defining 



w(r) 



m 



£+1 



w(x), any dependence of the eigenvalue equation (J3D on the mass m of the two 



bound-state constituents may be absorbed into a scaled (and therefore dimensionless) energy 
eigenvalue E and interaction potential V(x): 



E 



E 



V{x) 



V(r) 



m m 

The spinless Salpeter equation @ is then equivalent to the integro-differential equation || ||| 

oo 



E-V{x) x l+v w{x) 



7T 



dy G e (x,y)y 



£+1 







1 



d 2 2(^ + 1) d 
dy 2 



dy 



w{y) , (4) 



where the kernel Gi(x,y) is defined by 



G e (x, y) := 2* /+ 1 (~)' \ [(« " zf 2 K e (V^) - (s + zf 2 K t (v^ + i) 

with the abbreviations s = x 2 + y 2 and z = 2 x y. Here, Kg is the modified Bessel function of 
the second kind of order £ m . 



3 Conversion into an Equivalent Matrix Eigenvalue Problem 



In order to rephrase the eigenvalue equation (Q) in the form of an (easier-to-handle) algebraic 
problem, we expand every function on the positive real line R + we encounter into a complete 
orthonormal system {f n (x), n = 0, 1, 2, . . .} of basis functions for the Hilbert space L2{R + ). 
The solutions w(x) of Eq. (Q), in particular, are then obtained in the form 



w(x) 



N 

E 

n=0 



with some set of real coefficients c n . This treatment would be exact for N = oo and represents 
an approximation for N < oo of, however, increasing accuracy with increasing matrix size N. 
By application of this procedure, we are able to recast the spinless Salpeter equation (g) into 
the form of a matrix eigenvalue equation for the vector c = {c n } of the expansion coefficients: 



Ec 



(5) 
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Here, the elements of the "power matrices" and "potential matrices" are defined by 

oo 

P® := [dxx t+l f n (x)f m (x) = P£l, m,n = 0,1,..., N , 



o 

CO 



V® := / dx x e+1 V(x) f n (x) f m (x) = V%L, m,n = 0,1,..., N . 



o 

The matrices represent the action of the kinetic term on the vector of basis functions f n . 
The main advantage of our procedure is the fact that the "kinetic-energy matrices" have 
to be calculated only once (for a chosen matrix size N)\ The elements Tnm of these matrices, 
however, have to be calculated separately for every value of the orbital angular momentum I. 
The solution of the matrix equation (|5|) then gives the energy eigenvalues Ef., k = 0, 1, . . . , N, 
as well as the reduced radial wave functions 

TV 

» t (r) = m w ^ C f)/„(mr), k = 0,1, . . . ,N . (6) 

n=0 

Our choice of the basis functions f n (x) is primarily dictated by our demand to allow for an 
analytical treatment of the spinless Salpeter equation (|3|) to the utmost (reasonable) extent: 

f n (x) := \/2 exp(-x) L n (2x) , n = 0,1,2,... , 

where L n {x) are the Laguerre polynomials 

n / „ \ „t 

M*) = £("!)' 



n \ x" 
t ) tf 



i, - i 

t=o 

Trivially, these basis functions satisfy the orthonormalization condition 



d-X /n(S') fmipC) — $nm • 

For this choice of basis functions it is straightforward to write down the explicit expression for 
the matrix elements Pnm- 

This easy availability of the explicit expressions for the power matrices P^ strongly suggests 
to consider, as a special case, a class of interaction potentials of (generalized) power-law form: 

y(r) = £a n r", (7) 

where the two sets of (otherwise arbitrary) integers n and real constants a n (playing the role 
of coupling strengths) are only constrained by the necessary boundedness from below of the 
Hamiltonian H defined in Eq. ([!]): n > — 1 if a n < 0. For this class of interaction potentials, 
the potential matrices V^' simply become linear combinations of the power matrices P^ : 

y{i) = a n p (l+n) 

m n+l 

It goes without saying that the procedure described here works in the same way also, at least, 
for all power-law potentials involving additional exponential (damping) factors, that is, for all 
interaction potentials of the form 

V ( r ) = 2 a nr n exp(-6 n r) , b n > . 

n<=Z 
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4 The Kinetic-Energy Matrix 

Clearly, the main task in our game is the calculation of the kinetic-energy matrix T®. This 
enterprise is somewhat involved but, as already mentioned, has to be undertaken only once. 
In general, the elements of the kinetic-energy matrix are given by some expression of the 
form H 

2 n 

T nl = ~Y n, k) c{£; k,m) , m, n = 0, 1, . . . , N . 



£i=0 



Here, the factors D(£; n, k) represent the action of the differential operator on the right-hand 
side of the integro-differential equation (||) on the chosen basis functions f n (x). For our choice 
of basis functions, these factors read explicitly || 



D(t;n,k) = 2 V2^—!-2 k 
kl 



(k + £ + l) 



n 
k 



+ (k + 2£ + 2) 



n 



k + 1 



The factors c(£; k, m) are the expansion coefficients, in terms of the basis functions f n (x), of a 
particular integral over the kernel Gi(x,y) 0: 



c(£;k,m) = J dx f m (x) J dy G e (x,y) y e+k exp(-y) . 
o o 

For increasing orbital angular momentum £, the kernels Gp(x, y), £ = 0, 1, 2, . . . , become very 
rapidly rather complicated expressions. This circumstance forces us to calculate these factors 
c(£; k, m) separately for every single value of the orbital angular momentum £ = 0, 1, 2, ... we 
are interested in. Here, we confine ourselves to the discussion of the cases £ = and £ = 1. 
We find, for states with orbital angular momentum £ = (the so-called "S waves"), || 



c(0; k, m) 



[A/2] p m 

2V2k\ ]T ££(-i) 9+r 

p=0 q=0 r=Q 



k + 1 
2p+l 



P 

Q 



m 
r 



2 r (r + l){2k -2q + r + 1)! 



(2/c-2g + l)!!(2/c-2g + 2r + 3)!! ' 
and, for states with orbital angular momentum £ = 1 (the so-called "P waves") 

k + 1 



c(l; k, m) 



[A/2]+l p+l m 

2V2k\ x: EE(-!) 9+r 

p=0 q=0 r=0 



k + 2 
2p 



+ 



k + 1 
2p + l 



P 



1 



1 I (2k-2q + 5)l 



I P 
\ q 

m 
r 



1 (2A;-2g + 3)!! 

2 r (2k -2q + r + 4)! 
(2k-2q + 2r + 5)l\ 



where 



k 
2 

k-1 



for k even , 
for k odd , 



and 



(2n + 1)!! = 1 x 3 x ••• x (2n - 1) x (2n + 1) , n = 0,1,2, 



The calculation of the factors c(£; k, m) for £ > 1 is straightforward but somewhat involved 
because of the increasing complexity of the kernels Gg(x,y). 
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5 The Mathematica 4.0 Notebook Salpeter. 

The "semianalytical matrix method" developed here for the solution of the spinless Salpeter 
equation (||) has been implemented in a Mathematica 4.0 routine called Salpeter .nb (which 
may be obtained by contacting the authors). The routine Salpeter .nb requires as input, for 
two particles of mass m experiencing some interaction potential V(r) and for a chosen matrix 
size d = N+l, only this potential V{r) expressed in terms of the single terms a n r n in Eq. (|?]), 
these terms written in the form vpot [asubn , n , m , d] . For instance, for the "funnel" potential 

V(r) = + a\ r , 

r 

type 

v [m_,d_] : =vpot [-asubminusl , -1 ,m,d] +vpot [asubl , 1 ,m,d] . 

Upon entering the command etot [ell,m,d] , this routine then computes, for orbital angular 
momentum £, particle mass m, and matrix size d = iV+1, the bound-state energy eigenvalues 
Ek, k = 0, 1, . . . , N, together with the eigenvectors c^ k \ k = 0, 1, . . . , N, from Eq. (||) and the 
corresponding radial eigenfunctions Rk(?) = r Wk(r), k = 0, 1, . . . , N, according to Eq. (||). 

We illustrate the power of the "semianalytical matrix method" by applying it to the case 
of a harmonic-oscillator potential V(r) = ar 2 , a > 0, for the following reason: In momentum 
space, the operator r 2 is represented by the Laplacian w. r. t. the momentum p, r 2 — » — A p , 
while the kinetic energy T, nonlocal in configuration space, is represented by a multiplication 
operator. Only for a harmonic-oscillator potential the momentum-space representation of the 
semirelativistic Hamiltonian H is therefore of the form of a nonrelativistic Hamiltonian with 
a square-root interaction potential: 

H = -a A p + 2 yp 2 + m 2 . 

The resulting (nonrelativistic) Schrodinger equation is then solved with a standard numerical 
procedure designed exactly for this purpose || . (For details, in particular, for exact analytical 
upper and lower bounds on the energy levels, see, for instance, Refs. flQfl .) 



Table 1: Eigenvalues (in units of GeV) of the semirelativistic Hamiltonian H = 2 -\/p 2 + m 2 + V(r) 
with a harmonic-oscillator potential V(r) = ar 2 , for states of principal quantum number n — 1, 2, 3, 4 
and orbital angular momentum t = or £ = 1 (called nS or nP in the usual spectroscopic notation) 
obtained by the "semianalytical matrix method" for increasing matrix sizes d = 1, 3, 10, 25, compared 
with the corresponding "exact" eigenvalues obtained from the momentum-space representation of the 
spinless Salpeter equation, which, for the special case of a harmonic-oscillator potential, resembles a 
nonrelativistic Schrodinger equation (as discussed in the text). The chosen physical parameter values 
are m = 1 GeV for the particle mass and a = 0.5 GeV 3 for the harmonic-oscillator coupling strength. 



d = N + l 




£ 


= 






£ 


= 1 




IS 


2S 


3S 


4S 


IP 


2P 


3P 


4P 


1 


4.14531 








5.12166 








3 


3.91571 


7.08622 


11.94458 




5.03946 


7.68220 


14.80138 




10 


3.82522 


5.80930 


7.76609 


9.90246 


4.89944 


6.72710 


8.54739 


11.20168 


25 


3.82494 


5.79112 


7.48323 


9.01617 


4.90149 


6.69298 


8.28585 


9.74304 


exact 


3.82493 


5.79102 


7.48208 


9.00749 


4.90145 


6.69305 


8.28464 


9.74276 



Table [l] demonstrates the rapid convergence for increasing matrix size of the lowest-lying 
energy eigenvalues of the spinless Salpeter equation (pi) with a harmonic-oscillator potential. 
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